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I.  INTRODUCTION 


Improving  heat  dissipation  is  of  interest  to  many  applications,  especially  electronics.  As 
the  size  of  electrical  components  continues  to  decrease,  the  amount  of  heat  generated  in  the 
components  increase.  Composite  materials,  in  which  many  electronics  are  embedded,  have 
shown  poor  performance  dissipating  heat  generated  by  electrical  components.  For  instance, 
carbon  fiber  composites  have  gained  wide  use  in  automotive-to-munition-to-athletic  applications 
because  of  their  strength  and  light  weight;  however,  their  heat  dissipation  is  poorer  than  the 
metals  they  replace.  To  overcome  this  bottleneck,  techniques  like  Transient  Liquid  Phase 
Sintering  (TLPS)  and  nanoparticles  have  been  used  to  improve  heat  dissipation  from  carbon 
fibers.  To  improve  heat  dissipation  in  carbon  fiber  composites,  interactions  at  the  interface  of 
these  materials  require  further  understanding. 

Composites  combine  two  or  more  materials  with  different  properties  to  yield  a  material 
with  unique  properties  [1  ].  Nanocomposites  are  a  class  of  composites  which  contain  at  least  one 
material  with  at  least  one  nanoscale  dimension.  Nanoscale  fillers  have  properties  that  can  be 
different  from  their  bulk  material  due  to  size  effects.  In  studies,  nanocomposites  have  shown 
unproved  electrical,  mechanical,  and  thermal  properties.  Carbon  fiber  composites  have  shown 
success  in  improving  mechanical  properties.  Carbon  nano  tube  (CNT)-based  nanocomposites 
have  been  studied  for  integration  with  carbon  fiber  composites  to  improve  heat  transport. 

Carbon  fibers,  like  other  carbon  allotropes,  have  poor  thermal  transport  in  transverse 
directions  (out  of  and  across  fiber  planes)  but  excellent  transport  in  longitudinal  (along  the  fiber) 
directions.  In  most  applications,  carbon  fibers  are  combined  with  an  epoxy.  Epoxy  is  a 
thermo  setting  polymer  that,  like  the  carbon  fiber,  has  poor  thermal  conductivity;  however, 
carbon  fiber  reinforced  epoxy  composites  are  strong  and  lightweight.  Since  neither  carbon  fiber 
nor  epoxies  possess  high  thermal  transport,  several  studies  have  focused  on  how  to  improve  heat 
dissipation  in  the  carbon  fiber  composite.  Techniques  have  varied  based  on  the  method  of 
production  for  the  carbon  fiber  composites.  In  the  case  of  filament  wound  carbon  fibers 
composites,  efforts  have  focused  on  modifying  the  epoxy  which  acts  as  the  bonding  agent 
holding  the  fibers  together. 

Molecular  Dynamics  (MD)  simulations  can  be  used  to  track  the  trajectory  of  individual 
atoms  and/or  molecules.  By  applying  the  appropriate  interatomic  potential  to  govern  the 
interactions  of  the  atoms,  MD  simulations  can  be  used  to  understand  the  thermal  transport  at  the 
interface  of  materials.  Since  MD  simulations  calculate  atom  positions,  properties  such  as  kinetic 
energy,  dispersion,  and  thermal  conductivity  can  be  deduced. 

II.  THEORY 

Thermal  conductivity  can  be  calculated  using  MD  simulations  by  application  of  Fourier’s 
law  of  heat  conduction: 


1  AQ  dT 
A  At  dx 

where,  k  is  the  thermal  conductivity  of  the  system,  A  is  the  cross-sectional  area,  AQ/At  is  the  rate 
of  energy  transfer,  and  dT/dx  is  the  temperature  gradient.  The  area  and  simulation  time  of  the 
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system  are  known  values.  In Non-Equilibirum Molecular  Dynamics  (NEMD)  simulations,  two 
methods  can  be  used  to  generate  a  temperature  gradient  The  first  method  establishes  a  constant 
temperature  difference  in  the  system  and  calculates  the  energy  transferred  in  the  system. 
Alternatively,  the  energy  applied  to  the  system  can  be  constant  and  the  temperature  gradient 
calculated. 

The  constant  energy  method  is  used  in  this  report.  To  cause  a  temperature  difference,  a 
constant  amount  energy,  Q ,  is  added  and  subtracted  from  heat  baths.  The  addition  of  energy 
raises  the  temperature  of  the  slab.  Conversely,  subtracting  lowers  the  slab’s  temperature. 

Figure  1  shows  the  simulation’s  arrangement  where  the  heat  baths  are  located  at  0.25  and  0.75  of 
the  system’s  total  length  [2],  which  allows  a  temperature  gradient  to  be  established  without  the 
heat  baths  being  located  next  to  one  another.  The  energy  must  be  divided  by  half  for  calculating 
thermal  conductivity  because,  in  this  technique,  the  energy  can  flow  in  two  directions  from  the 
baths.  The  system  is  divided  into  slabs  of  equal  thickness  in  the  longitudinal  direction.  The 
temperature  gradient  is  calculated  based  on  statistical  mechanics.  If  the  kinetic  energy,  KE,  of  an 
atom  is  known,  then  the  temperature  of  that  atom  is  T=2(KE)/3ks,  where  ks  is  Boltzmann’s 
constant.  The  temperature  of  each  slab  is  calculated  by  averaging  together  the  temperatures  of 
the  atoms  within  that  slab.  For  Fourier’s  law  to  be  applicable,  the  temperature  gradient  produced 
must  be  linear. 


Figure  1.  NEMD  Setup  for  Constant  Energy  Temperature  Gradient 


III.  NANOCOMPOSITE  MODEL 

The  flowchart  of  a  MD  simulation  is  shown  in  Figure  2.  The  initial  step  in  the  simulation 
is  to  assign  the  coordinates  for  each  atom  Two  software  packages  were  used  to  generate  the 
initial  coordinates  for  the  nanocomposite:  TubeGen  and  Packmol  [3,  4].  TubeGen  was  used  to 
generate  the  initial  configuration  for  CNTs.  An  example  TubeGen  script  is  shown  in  the 
appendix  which  would  produce  the  coordinates  for  a  (10,  10)  Single-Wall  Nanotube  (SWNT), 
which  is  frequently  used  in  MD  simulations  [5].  Figure  3  shows  a  cartoon  representation  of  the 
coordinates  using  Visual  Molecular  Dynamics  (VMD)  [6]. 
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Figure  2.  Description  of  MD  Simulation  Algorithm  [7] 


Figure  3.  VMD  Rendering  of  (10,10)  SWNT From  Coordinates  Generated  by  TubeGen 

Similarly,  Packmol  can  be  used  to  generate  the  polymer  molecules’  initial  coordinates  by 
taking  the  coordinates  of  a  typical  molecule,  like  ethylene  (Figure  4),  and  finding  the  most 
efficient  packing.  Using  the  Protein  Data  Bank  (PDB)  [8],  the  coordinates  of  myriad  atoms  can 
be  found.  An  example  script  for  Packmol  to  pack  ethylene  molecules  into  a  80-by-80-by- 100 
nanometer  box  with  a  cylindrical  void  in  the  center  to  house  the  CNT  is  shown  in  the  appendix. 
Since  a  Large-Scale  Atomic  and  Molecular  Massively  Parallel  Simulator  (LAMMPS)  is  used  as 
the  MD  simulator  [9],  the  coordinates  must  be  formatted  for  use  in  LAMMPSs.  VMD  has  a  set 
of  tools  (TopoTools)  that  can  be  used  to  generate  a  LAMMPS-readable  format  [6], 
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Figure  4.  Ethylene  Monomer  Produced  From  Coordinates  in  PDB 
and  Rendered  Using  VMD  (Green,  C  Atoms;  Gray,  H Atoms) 

To  achieve  meaningful  results  from  a  MD  simulation,  the  appropriate  interactions  must 
occur  between  the  atoms.  The  forces  between  atoms  are  governed  by  the  interatomic  potential. 
Several  interatomic  potentials,  also  known  as  force  fields,  exist.  The  differences  in  potentials 
arise  from  how  the  parameters  of  the  potential  were  fitted.  Some  potentials  are  based  on 
experimental  data,  whereas  others  come  from  quantum  mechanics  calculations.  The  Tersoff  and 
Reactive  Empirical  Bond  Order  (REBO)/ Adaptive  Inter  molecular  Reactive  Empirical  Bond 
Order  (AIREBO)  potentials  are  used  regularly  for  MD  simulations  involving  thermal 
conductivity  and  CNTs  since  they  describe  covalent  bonds  [10].  For  the  present  simulations, 
AIREBO  will  be  used  since  it  addresses  the  interactions  of  not  only  carbon  (C)  atoms,  but  also 
hydrogen  (H)  atoms.  Additionally,  the  van  der  Waals  interaction  of  the  CNT  with  the  polymer 
and  ofpolymer  chains  is  described  by  the  Lennard-Jones  potential. 

Finally,  simulations  are  performed  by  advancing  the  atoms  and  recalculating  the  forces.  A 
time  step  is  chosen  that  provides  stable  simulations,  yet  does  not  require  a  large  amount  of 
computational  tune.  Time  steps  on  the  order  of  a  femtosecond  are  usually  sufficient  for  tracking 
phonons;  however,  the  H  atoms  can  sometimes  cause  high  frequencies  that  require  smaller  tune 
steps.  Simulations  usually  require  an  equilibration  period  to  allow  the  total  energy  of  the  system 
to  stabilize.  After  stabilization,  the  desired  properties  can  be  calculated. 

IV.  SIMULATION  SETUP 

The  (10,10)  SWNT  and  polyethylene  (PE)  are  used  to  model  nanocomposite  for  thermal 
conductivity  measurements.  The  (10,10)  SWNT  is  chosen  because  the  diameter  of 
1.4  nanometers  is  reported  as  a  typical  size  for  SWNTs  [5].  PE  is  a  simple  hydrocarbon,  yet  it  is 
applicable  since  polymers  generally  have  low  thermal  conductivities.  Its  purpose  is  to  show  the 
trend  of  interacting  with  the  CNT.  The  CNT  and  PE  are  generated  for  a  length  of 
100  nanometers. 
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The  Lennard- Jones  potential  describes  the  interactions  of  non-bonded  atoms.  The 
Lennard- Jones  potential,  U,  is  given  by  the  following: 


U  =  4s 
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where,  e  is  the  depth  of  the  potential  well,  a  is  the  distance  at  which  the  potential  energy  is  zero, 
and  r  is  the  distance  between  two  atoms.  The  parameters  a  and  s  depend  on  the  types  of  atoms 
being  considered.  For  CNTs,  the  Lennard- Jones  parameters  are  the  same  as  for  graphite  (cr  =  3.4 
A  and  e  =  4  meV).  The  Lennard- Jones  parameters  for  PE  are  based  on  Reference  1 1  where  the 
values  for  C  and  H  atoms  are  listed.  When  the  two  atoms  are  the  same,  then  Lennard- Jones 
parameter  values  can  be  used  as  listed;  however,  for  atoms  that  are  different  (or  in  the  case  ofC 
atoms  that  have  different  hybridizations),  the  interaction  must  be  calculated  using  the  Lorentz- 
Berthelot  mixing  rules: 

£ij  ~  yf^ii^jj 
and 

aij  =2^°u  + 


where,  i  and  j  are  the  atom  subscripts. 

Simulations  are  performed  using  LAMMPS  simulation  software.  Periodic  boundary 
conditions  are  maintained  in  all  directions.  The  time  step  used  for  minimization  is 
0.1  femtosecond.  Each  system’s  energy  is  first  minimized  using  the  conjugate  gradient  method. 
The  appendix  contains  a  sample  of  the  minimization  script. 

V.  RESULTS 

As  a  comparison  to  Reference  12,  a  PE  matrix  is  modeled  to  begin  the  proof  of  concept. 
Minimization  is  performed  on  the  polymer  matrix  using  the  minimization  command  available  in 
LAMMPS.  Figure  5  shows  the  effectiveness  of  energy  minimization  to  make  a  more  compact 
structure.  Green  represents  the  C  atoms,  and  gray  represents  the  H  atoms.  Minimization,  which 
finds  the  configuration  with  the  lowest  potential  energy  for  the  system,  is  a  necessary  step  since  a 
simulation  that  starts  at  a  high  energy  will  likely  become  unstable  and  the  system  will  “blow  up.” 
PE  in  a  volume  of  40  nanometers  squared  by  100  nanometers  in  length  is  chosen.  The 
minimized  structure  is  shown  in  Figure  6,  along  with  a  plot  of  its  energy  minimization  versus 
time  step.  The  energy  begins  as  a  positive  value  since  the  forces  are  calculated  between  each 
pair  of  atoms  using  the  interatomic  potential.  The  minimization  algorithm  moves  the  atoms  a 
small  time  step  to  incrementally  adjust  the  position  of  the  atoms  to  an  equilibrium  position, 
which  is  not  too  close  to  cause  repulsion  nor  too  far  to  cause  attraction.  The  largest  reductions  in 
energy  occur  in  the  first  few  time  steps.  Similarly,  the  before  and  after  minimization  structures 
of  the  composite  are  shown  in  Figure  7.  The  composite’s  volume  is  80  nanometers  squared  by 
100  nanometers,  which  allows  the  inclusion  of  the  SWNT.  Since  the  CNT  fills  the  center  of  the 
PE  matrix,  the  compact  minimization  does  not  occur  like  with  the  PE  alone.  In  the  case  of  the 
composite,  the  minimization  efforts  are  less  obvious;  however,  the  decreasing  total  energy  in 
Figure  8  verifies  that  minimization  did  occur. 
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(a)  Before  Energy  Minimization  (b)  After  Energy  Minimization 

Figure  5.  Minimization  Results  Shown  for  80  Ethylene  Molecules 
(Green,  C  Atoms;  Gray,  H Atoms) 


Figure  6.  Compacted  Block  (40-by-40-by-100  Nanometers) 
of  PE  After  Total  Energy  Minimization 


(a)  Before  Minimization  (b)  After  Minimization 

Figure  7.  PE/SWNT  Composite 
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Figure  8.  Energy  Minimization  of  Total  Energy  of  PE/SWNT  Composite 

VI.  CONCLUSIONS 

Several  strides  were  made  in  studying  the  PE  and  nanocomposite  models.  An  initial 
configuration  was  able  to  be  established  for  PE  and  a  PE/SWNT  composite.  Since  an 
appropriate  interatomic  potential  was  chosen  (the  AIREBO  potential),  minimization  was 
performed  on  both  systems.  Minimization  of  the  PE  structure  did  result  in  a  more  compact  lower 
energy  configuration;  however,  both  structures’  total  energy  decreased  significantly  during  the 
minimization  process.  The  modeling  of  polymers  is  a  computationally  expensive  process  which 
requires  more  computing  power  than  the  office  desktop  system  that  was  available  for  the  project. 

Ml.  FUTURE  DIRECTIONS 

Due  to  tune  constraints,  additional  tasks  were  not  able  to  be  completed  but  nonetheless  are 
needed  to  understand  thermal  transport  at  the  composite  interface.  Additional  models  that  need 
to  be  considered  include  arrays  of  SWNTs  since  the  actual  thermal  interface  incorporates  aligned 
CNTs.  To  extend  this  concept  further,  Multi-Walled  Carbon  Nano  tubes  (MWNTs)  should  be 
modeled  as  well  and  in  an  arrayed  geometry.  The  typical  CNT  is  multi- walled;  therefore, 
double-  or  triple-walled  CNTs  should  be  useful  in  establishing  trends.  The  actual  thermal 
conductivity  measurements  were  performed  using  the  constant  energy  NEMD  method  mentioned 
in  Section  II. 


80  nm  x  80  nm  x  100  nm  SWNT/PE 


7/8  (Blank) 
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A I  RE BO 

LIST  OF  ABBREVIATIONS,  ACRONYMS,  AND  SYMBOLS 

Adaptive  Inter  molecular  Reactive  Empirical  Bond  Order 

C 

carbon 

CNT 

carbon  nanotube 

H 

hydrogen 

LAMMPS 

Large-Scale  Atomic  and  Molecular  Massively  Parallel  Simulator 

MD 

Molecular  Dynamics 

MWNT 

Multi-Walled  Carbon  Nanotube 

NEMD 

Non-Equalibrium  Molecular  Dynamics 

PDB 

Protein  Data  Bank 

PE 

polyethylene 

REBO 

Reactive  Empirical  Bond  Order 

SWNT 

S  ingle- Wa  11  Nanotube 

TLPS 

Transient  Liquid  Phase  Sintering 

VMD 

Visual  Molecular  Dynamics 
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APPENDIX 
EXAMPLE  SCRIPTS 


A.l  Example  TubeGen  Script 

set  relax_tu be  yes 
set  format  xyz 
set  units  angstrom 
set  elementl  C 
set  element2  C 
set  bond  1.4210 
set  shape  hexagonal 
set  chirality  10,10 
set  gutter  10, 10,  0 
set  cell_count  1, 1,  41 
generate 
save  cnt.xyz 
exit 


A.2  Example  Packmol  Script 

tolerance  2.0 
seed  1211984 

output  polyethylen  e_cylin  drica  l_vo  id.  p  db 
filetype  pdb 
nloop  1000 

structure  ethylene_10.pdb 
number  266 

inside  box  0.  0.  0.  80.  80.  100. 
outside  cylinder  40.  40.  0.  0.  0. 1.  15.  100. 
end  structure 


A.3  Example  LAMMPS  minimization  Script 

clear 

log  log.polyethylene_minimize 
### Minimization  of  PE#### 
units  metal 

boundary  p  p  p 

atom_style  atomic 

pair_style  airebo  3.011 

read_data  crys  ta  I/p  olyethylen  e_box 

pair_coeff  *  *  ../LAMMPS/lammps-171unll/potentials/CH.airebo  C  H 

neighbor  2.0  bin 
neigh_modify  every  3  delay  3 
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time  step  0.0001 


#####minimiza  von##### 

dump  1  all  custom  5000  /scratch/walkerel/dump. polyethylene  id  typex  y  z 
dump  2  aiixyz  5000 /scratch/walkerel/dump. polyethylene. xyz 
thermo  2 

thermo_style  custom  step  temp  ke  etotal pe  vol 
minimize  0.01.0e-8  1  00000  1  00000 
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